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ABSTRACT 

Using our recently improved Monte Carlo evolution code, we study the evolution of the binary 
fraction in globular clusters. In agreement with previous ./V-body simulations, we find generally that 
the hard binary fraction in the core tends to increase with time over a range of initial cluster central 
densities for initial binary fractions < 90%. The dominant processes driving the evolution of the core 
binary fraction arc mass segregation of binaries into the cluster core and preferential destruction of 
binaries there. On a global scale, these effects and the preferential tidal stripping of single stars tend to 
roughly balance, leading to overall cluster binary fractions that are roughly constant with time. Our 
findings suggest that the current hard binary fraction near the half-mass radius is a good indicator 
of the hard primordial binary fraction. However, the relationship between the true binary fraction 
and the fraction of main-sequence stars in binaries (which is typically what observers measure) is 
non-linear and rather complicated. We also consider the importance of soft binaries, which not only 
modify the evolution of the binary fraction, but can drastically change the evolution of the cluster as 
a whole. Finally, we describe in some detail the recent addition of single and binary stellar evolution 
to our cluster evolution code. 

Subject headings: globular clusters: general — methods: numerical — stellar dynamics 



1. THE BINARY FRACTION 

Observations and recent theory strongly suggest that 
the initial mass function ( IMF) is univers al among non- 
zero me tallicity stars (e.g. . IChabrierl 12003). Indeed. If3atd 
(|2008bD suggests that radiative feedback may naturally 
regulate the star formation process so as to produce 
an IMF that is only weakly dependent on the prop- 
erties of the progenitor molecular cloud. Naively, one 
would also expect that other features of the initial stellar 
population — like the binary fraction — should be nearly 
universal. Hydrodynamical star formation simulations 
yield companion star frequencies and binary fractions 
that are largely independent of the properties of the pro- 
genitor molecular cloud (although the statistics in some 
cases are marginal), and are quite consistent with ob- 
servations (|Batdl2008aHBate et al.ll2003tlBate fe Bonnelll 
HI). 

Observations of stars in low stellar density environ- 
ments where dynamics is unimportant, such as the solar 
neighborhood, yield a binary fraction of ~ 50% among 
solar-type stars, with an increasing trend with primary 
mass (e.g.. iDuauennov fe Mavorill991l : iFischer fc Marcvl 
Il992t) . O pen clusters sim ilarly show such large binary 
fractions (|Fan et al.l 1 1996f ) . However, observations of 
dense globular cluster cores typically yield binary frac- 
tions that are significantly smaller. HST observations of 
the core-collapse cluster NGC 6397 yield a binary frac- 
tion of ~ 5% in the core a nd « 1% beyond the half-mass 
radius l|Davis et al.l l2008h . For the canonical non core- 
collapse cluster 47 Tuc, the binary fraction is f=a 13% 
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(Alb row eEaDISOQil)- The core binary fraction generally 
ranges from a few percent to tens of percen t , approaching 
50% in some cases for less dense clusters (|Sollima et al.l 
l2007f ). Where measured, the binary fraction outside th e 
core is always smaller (see table in iDavis et al.l l2008| ) . 
The question naturally arises: Are the currently observed 
relatively small core binary fractions in globular clusters 
consistent with initially larger binary fractions of ~ 50%? 

There are many strongly coupled processes that deter- 
mine the evolution of the core binary fraction in a dense 
stellar system. Stellar evolutionary processes alone can 
affect the properties of a binary greatly, causing it to ex- 
pand or shrink via mass transfer or winds, circularize via 
dissipative effects, lose mass, receive a systemic velocity 
kick due to a supernova, or disrupt or merge. The prop- 
erties of the binary feed into the dynamical interaction 
rate with other stars or binaries, causing it to interact 
more or less frequently depending on its semimajor axis, 
eccentricity, mass, and systemic velocity. A strong dy- 
namical interaction of a binary can disrupt it, exchange 
one of its members for an incoming star, cause its orbit 
to expand or shrink, modify its eccentricity, increase its 
systemic velocity via gravitational recoil, or cause two or 
more stars to physically collide. The dynamically mod- 
ified binary properties feed back into binary stellar evo- 
lution, possibly initiating or halting mass transfer, or in- 
creasing tidal effects. In contrast to stellar evolutionary 
processes, the dynamical interaction rate depends on the 
cluster density and velocity dispersion, which evolve with 
time. Since binaries are typically more massive than sin- 
gle stars, mass segregation can increase their numbers in 
the core at the expense of single stars. The tidal effects of 
the host galaxy will preferentially strip single stars from 
the halo of the cluster. 

For a globular cluster of typical mass (~ 10 5 Mq) and 
size (half- mass radius ~ 3pc), its global evolution can 
be divided into three phases according to the timescales 
of the relevant physical processes. At early times (~ 
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few x 10 Myr), the evolution is largely driven by stellar 
evolutionary mass loss from the most massive stars in 
the cluster. At intermediate times (~ few Gyr), as mass 
loss from stellar evolution has slowed, the evolution is 
driven primarily by two-body r elaxation. At late times 
(possibly beyond a Hubble time; [Hurley 2007), when the 
core has reached sufficiently high density for binaries to 
strongly interact dynamically and release enough energy 
to prevent core collapse, the properties of the cluster are 
determined by the makeup of the binary population in 
this quasi-equilibrium "binary burning" phase. 

The core binary fraction is clearly a quantity that is 
affected by nearly all physical processes operating in a 
cluster, and is of obvious observational interest. Compar- 
ing observed core binary fractions with simulation results 
(in combination with other obscrvables) is thus a good 
measure of our theoretical understanding of cluster evo- 
lution. There can be dramatic differences in definition 
between the observed binary fraction and what theorists 
call the binary fraction, however. 

When measured with the common offset main- 
sequence (MS) method, MS-MS binaries are detected by 
their appearance as distinctly brighter MS objects. The 
observed binary fraction is defined as the ratio of the 
number of these "binary sequence" objects to the total 
number of objects in the MS and the binary sequence, 
corrected for the assumed number of binaries with mass 
ratio so small they would blend in with the MS. 

The theorists' definition of the binary fraction is typi- 
cally the ratio of the number of binaries to the total num- 
ber of "objects" (single stars or binaries). Furthermore, 
computational theorists tend to consider only "hard" bi- 
naries. That is, binaries with binding energy greater 
than the typical particle energy, which typically become 
more tightly bound (harden) as a result of encounters 
(jHeggie fc Hudl20M ). Soft binaries— b inaries with bind- 
ing energy less than the typical particle energy in a clus- 
ter, which typically become less tightly bound (soften) 
or dissociate completely — are less frequently considered. 
We consider in detail the difference between the obser- 
vational and theoretical definitions of the binary fraction 
below, as well as the importance of soft binaries. 

Recently, two very different simulation methods have 
been used to study the e volution of the binary frac- 
tion. Ilvanova et alj (|2005T ) have developed a simplified 
Monte Carlo method in which a dense, massive cluster 
is modeled as a constant-density core plus halo (to sim- 
ulate the long-lived binary burning phase that clusters 
may reach late in their evolution). Binaries and stars 
are evolved via the pop ulation synthesis code Star Track 
([Belczynski et al. 2008), and the strong dynamical inter- 
action s of binaries are inte grated numerically with Few- 
body (|Freeeau et all 120041 ). Objects move between the 
core and the halo due to mass segregation and systemic 
velocity changes resulting from dynamical encounters. In 
this approach the core mass increases slowly with time, 
with very few stars leaving the core after mass segregat- 
i ng into it. 

Ilvanova et al.1 (|2005| ) find, generally, that the core bi- 
nary fraction decreases significantly with time. Even for 
a modest core density of 10 3 pc -3 , they find that an ini- 
tial binary fraction of 100% yields a core binary fraction 
of 27% at 14 Gyr. For the density of 47 Tuc, they find 
that a 100% initial binary fraction yields an 8% core 




Fig. 1.— Evolution of iV-bodv (e.g.. IHurlev et al J I2007T > and 
simplified Monte Carlo lllvanova e t al. 2005) cluster models in core 
number density — binary fraction space. Each model's evolution is 
represented as a simple arrow, with the tip at the final properties, 
and the tail at the initial properties. For simplified Monte Carlo 
the final properties are measured at an age of f4 Gyr. For ./V-body 
they are measured at ~ 15 Gyr in most cases, with the 50% initial 
binary fraction model measured at 4 Gyr, and the 10 3,5 pc~ 3 core 
density model measured at 9 Gyr. Note that the binary fractions 
plotted here include only hard binaries. For reference, we plot as 
open circles the current observed properties for several Galactic 
globular clusters where measuremen t is possible, with data taken 
from the table in lDavis et al.1 1120081 1. 

binary fraction at 14 Gyr. It should be noted, how- 
ever, that these figures include substantial numbers of 
soft binaries — binaries that are so wide they are quickly 
destroyed by dynamical encounters. If only the hard bi- 
naries in these simulations are counted, an initial binary 
fraction of 25% in a 10 3 pc -3 core density cluster yields 
a 15% core binary fraction at 14 Gyr. For a density of 
10 5 pc -3 the core binary fraction evolves from an initial 
2 5% to 7%. 

IHurlev et all (|2007f ) have used a direct TV-body 
method, coupled wi th the BSE single and bin ary stellar 
evolution routines (|Hurlev et al.l l2000l l2002f ), to study 
the evolution of the binary fraction. The great benefit of 
this method is that it makes no simplifying assumptions 
about the underlying cluster evolution. On the other 
hand, it is computationally expensive, currently limiting 
its application to ope n clusters or g l obular s with low ini- 
tial binary fractions. IHurlev et al.l (|2007f ) find that the 
core (hard) binary fraction generally increases with time. 
For a cluster of 5 x 10 4 stars with a central density of 
~ 10 3 5 pc -3 , the core binary fraction rises from an ini- 
tial 20% to 52% at 9 Gyr. For lower initial densities the 
degree of increase of the core binary fraction is similar. 

On the face of it, the discrepancy between the two 
methods appears irreconcilable. However, the two meth- 
ods operate at very different core densities and cluster 
masses, both of which affect the half-mass relaxation 
time and hence the mass segregation timescale, as well as 
the binary dynamical interaction rate. Fig. [T] shows the 
evolution of the various models in core number density — 
binary fraction space. Note that the binary fractions 
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Fig. 2. — Evolution of the core to half-mass radius ratio for 
TV = 10 J initial models with 0%, 5%, and 10% primordial b inaries, 
comp aring our new MC results with those of direct TV-body jHurlevI 
120071 ). The thick solid lines show the TV-body models, with color 
denoting initial binary fraction, /;,. The thin solid lines show our 
MC simulation with all relevant physics turned on (stellar evolu- 
tion in singles and binaries, physical collisions, binary interactions, 
and a tidal boundary), again with color denoting the initial binary 
fraction. For the sake of comparison, MC simulations with stellar 
evolution turned off, and without a tidal boundary are shown in 
the thin dashed and dotted lines, respectively. For clarity, only 
the 0% initial binary fraction runs for these comparison models 
are shown, since the 5% and 10% do not differ appreciably from 
the 0% case. Clearly the evolution of this model is driven primar- 
ily by the effects of stellar evolution. With the exception of the 
increased expansion of the cluster core at early times in the MC 
model, there is very good agreement between MC and TV-body for 
all three binary fractions considered, suggesting that the imple- 
mentation of stellar evolution in the MC code is consistent with 
that of TV-body. Note that for the sake of comparison, the core 
radius here i s calculated using the s tandard definition for TV-body 
simulations (Cascrtano & Hut 1985). 

plotted here include only hard binaries. Each model's 
evolution is represented as a simple arrow, with the tip at 
the final properties, and the tail at the initial properties. 
It is clear from this figure that the two methods repre- 
sent very different regions of parameter space, and could 
simply be displaying different aspects of the same under- 
lying physics. The only point of concern is the TV-body 
model starting at ~ 10 pc~ 3 and evolving toward a 
higher binary fraction, nestled between two Monte Carlo 
models evolving in the opposite direction. 

To elucidate the evolution of the binary fraction, and 
to address the discrepancy between the existing TV-body 
and simplified Monte Carlo models, we have performed 
a grid of simulations with our newly-upgraded Monte 
Carlo cluster evolution code. Note that our Monte Carl o 
code is very different from that of llvanova et al.l (|2005l ). 
While their code assumes a constant core density with 
time and samples binary interactions using Monte Carlo 
techniques, our code sclf-consistcntly models the global 
evolution of a cluster, using Monte Carlo techniques to 
sample the stellar distribution function when applying 
the effects of two-body relaxation. The naming clash is 
the unfortunate consequence of the popularity and ap- 
plicability of Monte Carlo techniques in general. 

2. MODERN SIMULATIONS 

Our Monte Carlo (MC) code self-consistently models 
the evolution of star clusters due to the effects of two- 
body relaxation, evaporation through a Galactic tidal 
boundary, dynamical scattering interactions of binaries, 



physical stellar collisions, and now single and binary stel- 
lar evolution. The details of the method an d its imple- 
ment ation are described in detail elsewhere (Ijoshi et ail 
120001 f200ll : iFreeeau et aD 120031 : iFreeeau fc Rasioll2007h . 
Here we focus on the addition of stellar evolution. 

For coding simplicity and for more directed compar- 
isons with existing TV-body simulations, we have incor- 
porated the BSE single and bi nary stellar evolution rou - 
tincs in our Monte Carlo code (jHurlev et alJl2000L[2002T l. 
In our code stellar evolution is performed for each object 
(single star or binary) during a timestep in step with 
dynamics. Since at early times a cluster can lose a lot 
of mass due to supernovae, we make sure to limit the 
timestep so that no more than a small fraction of the to- 
tal cluster mass is lost in one step (typically we set this 
fraction to 10~ 3 ). 

To test that our inclusion of the stellar evolution rou- 
tines is accurate, we have compared with the V-body re- 
sults of lHurlevI (|2007f ). who evolved V = 10 5 cluster mod- 
els with binary fractions ranging from to 10%. The re- 
sults arc shown in Fig. which displays the evolution of 
the core to half-mass rad ius ratio (r c /rh) with time. The 
data from iHurlevi (|2007t l were extracted from that pa - 
per using ADS's Dexter applet (|Demleitner et al.ll200i[) . 
For reference, we also plot the evolution of a model with 
stellar evolution turned off, and a model without an ex- 
ternal tidal field. Since the model without stellar evo- 
lution reaches core collapse in under 1 Gyr, and since 
the model with no tide differs only minimally from the 
models with all physics turned on, it's clear that stellar 
evolution drives the evolution of the cluster. These mod- 
els are thus a good test of our treatment of stellar evo- 
lution. The agreement with V-body is quite good given 
the vastly different methods, although the Monte Carlo 
models tend to expand more at early times due to super- 
novae. The peculiar feature that the evolution doesn't 
appear to depend strongly on initial binary fraction is 
reproduced in our models. At late times (> 15 Gyr) our 
models begin to diverge with V-body. This is likely due 
to the fact that the clusters have lost roughly 70% of their 
stars by this time, and our apocenter-based treatment of 
the tide will tend to underestimate tidal mass loss as the 
number of cluster stars decreases, when an energy-ba sed 
criterion is more appropriate fe.g.. iGiersz et a l. 2008). 



In a forthcoming paper we will perform more detailed 
comparisons with existing models of the open cluster 
M 67 and the g lobular clusters M 4 and NGC 6397 



in the literature (|Hur 


ev et all 120051: IGiersz et all 120081: 


iHeggie & Giersz 2008: 


Giersz & Heeeie 20091). Given the 



vast differences between the V-body method and our 
Monte Carlo method, we take the a greement between 
our models and those of lHurlevI l)2007f ) in Fig. [5] as a sign 
that our implementation of BSE in our code is at least 
consistent with that in V-body. 



There is one aspect of our method that deserves special 
mention, however. It is generally believed that if a clus- 
ter avo ids a collisional runaway phase (e.g.. lFreitag et aTl 
I2006bllah the stellar -mass black holes formed early in a 
cluster's lifetime will quickly sink to the core and dynam- 
ically decouple from the rest of the cluster, undergoing 
their o wn evolution, much like an inde pendent small star 
cluster (|Sigurdsson fc Hernquistl 119931 ). The BH subsys- 
tem will quickly dissolve through its own internal dynam- 
ics, ejecting all but one or two of the BHs on a timcscalc 
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Fig. 3. — Evolution of our Monte Carlo cluster models in core 
number density — binary fraction space. Conventions are as in 
Fig. [T] Solid arrowheads represent values measured at 14 Gyr, 
while dotted arrowheads are values measured before tidal disrup- 
tion (since these models didn't last for 14 Gyr) at times between 
~ 8 and ~ 13 Gyr for the medium initial density models, and be- 
tween ~ 2 and ~ 12 Gyr for the high initial density models. For 
reference, the detailed evolution of the low-density f b = 0.05 model 
is shown in small gray dots. The low initial density models have 
initial half-mass relaxation times of = 0.8 Gyr, the medium 
density models have t r h = 0.3 Gyr, and the high density models 
have t rh = 0.09 Gyr. 



< 1 Gyr. Aside from removing nearly all BHs from the 
cluster, the result is a mild energy injection into the clus- 
ter core, causing it to expand somewhat at early times 
(|Mackev et al.ll2007t ). A typical star cluster of N = 10 6 
objects will contain a subsystem of up to re 10~ 3 iV = 10 3 
BHs evolving independently in the core (jO'Learv et alJ 
l2006f ). For the N = 10 5 clusters considered in this work, 
the number is re 100. The Monte Carlo method is not 
designed to handle subsystems of less than a few hun- 
dred objects, since they are often far from spherically 
symmetric, and large angle scattering dominates. (Note 
that we treat small- N encounters up to N — 4 via direct 
integration.) We therefore truncate the mass function at 
18.5 Mq (the largest progenitor mass not resulting in a 
BH) for the runs presented here. The resulting discrep- 
ancy in r c /rh is important only at early times and as 
Fig. [5] shows is minimal. 

We have performed several simulations of evolving 
clusters for a grid in initial binary fraction and initial 
cluster virial radius (or equivalently, central density). All 
our simulations start with N = 10 5 objects initially (an 
object being ei ther a binary o r a sin gle star), and like the 
simulations of iHurlev et alJ (|2007l ). assume a Plummer 
density profile with no primordial mass segregation, a 
"standard" Galactic tide (cluster at 8.5 kpc from Galactic 
center , 10 11 Mq Galactic mass enclosed), a lKroupa et alJ 
l| 1993ft IMF, and only hard binaries. Our IMF extends 
from 0.15 to 18.5 Mq, binary secondary masses are drawn 
from a distribution flat in the mass ratio, the semimajor 
axis a is drawn from a distribution flat in logo from a 
minimum of a m i n = 5(i?i + R2), where Ri are the indi- 
vidual stellar radii, to a maximum corresponding to an 
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Fig. 4. — Evolution of our Monte Carlo cluster model starting 
with = 0.05 and n c £S 10 2,5 pc -3 . The top panel shows the 
evolution of r c /r h with time. The cluster enters a binary burning 
phase at ~ 12 Gyr. The bottom panel shows the evolution of the 
half- mass radius of single stars, rh, s (solid line), and the half-mass 
radius of binaries, t, (dashed line), relative to the overall cluster 
half-mass radius, r^. The differential mass segregation between the 
single and binary populations is evident, with the single stars ex- 
panding slightly relative to the bulk of the cluster, and the binaries 
contracting significantly. The quantity b/ r h decreases steadily 
until ~ 11 Gyr due to mass segregation, at which point it begins to 
increase due to destruction of binaries preferentially in the cluster 



orbital velocity of the lighter member equal to the local 
velocity dispersion, and the eccentricity is drawn from a 
thermal distribution truncated at the value correspond- 
ing to contact at a m ; n . Note that our large a cutoff for 
wide binaries is e quivalent to the hard -soft boundary for 
equal- mass stars (|Fregeau et al.|[2b06T ). 

Fig. [3] shows the evolution of our models in core num- 
ber density — binary fraction space. It is clear that for 
all the but the highest initial binary fraction cases, the 
core binary fraction increases with time. The observa- 
tional data points seem to be consistent only with cluster 
models that started with relatively low central densities 
(~ 10 3 pc~ 3 ) and small hard binary fractions (~ 5%). 
As we discuss in the next section, the core binary frac- 
tion is typically estimated observationally by measuring 
the fraction of main sequence stars belonging to the bi- 
nary main sequence, and convolving it with an assumed 
binary mass ratio distribution. It is not a priori evident 
that this MS binary fraction reflects the underlying true 
binary fraction. 

Why does the core binary fraction generally tend to in- 
crease with time? As mentioned above, there are many 
strongly coupled processes that affect the core binary 
fraction. However, the general trend can be understood 
approximately as an interaction between mass segrega- 
tion of binaries into the core, and the destruction of bi- 
naries preferentially in the core. 

Fig. [4] shows the evolution of our Monte Carlo clus- 
ter evolution model starting with ft, = 0.05 and n c re 
10 25 pc~ 3 . As the evolution of r c /rh in the top panel 
shows, the cluster core contracts steadily until it enters 
a phase of binary burning at the relatively late time of 
~ 12 Gyr. The bottom panel shows the evolution of the 
half-mass radius of single stars, r^ s (solid line), and the 
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Fig. 5. — Evolution of the number of single stars in the core, 
7V CjS (solid line), and number of binaries in the core, N c t, (dashed 
line), in our Monte Carlo cluster evolution model starting with 
ft = 0.05 and n c si 10 2 ' 5 pc -3 . The quantity N c ^ s declines steadily 
due to standard gravothermal evolution, in which the cluster core 
becomes denser and smaller in number with time. The quantity 
N c b is roughly steady until ~ 11 Gyr due to mass segregation of 
binaries into the cluster core. 

half-mass radius of binaries, r^b (dashed line), relative 
to the overall cluster half-mass radius, r/j. The differ- 
ential mass segregation between the single and binary 
populations is evident, with the single stars expanding 
slightly relative to the bulk of the cluster, and the bina- 
ries contracting significantly. The quantity rh.b/fft, de- 
creases steadily until ~ 11 Gyr due to mass segregation. 
It then begins to increase due to destruction of bina- 
ries preferentially in the cluster core by strong dynami- 
cal interactions and perturbed stellar evolution (see e.g., 
Ilvanova et alj 12005, for a discussion of perturbed binary 
evolution) . 

Fig. [5] shows the evolution of the number of single stars 
in the core, iV CjS (solid line), and number of binaries in the 
core, iV Cj b (dashed line) for the same model. The quan- 
tity N c , s declines steadily due to standard gravothermal 
evolution, in which the cluster core becomes denser and 
small er in number with time (e.g., iBinnev fc Tremainel 
2008). The quantity iV Cl b, on the other hand, is roughly 
steady until ~ 11 Gyr due to mass segregation of binaries 
into the cluster core. 

As suggested by Fig. [51 the core mass decreases with 
time, as expected from standard gravothermal evolu- 
tion. This is in contrast with th e simplified Monte Carlo 
method of Ilvanova et alj (|2005f ). in which the core mass 
steadily increases in time, due primarily to mass segre- 
gation of binaries into a core of fixed density. 

We note also that mass segregation of a binary into 
the core implies, by energy conservation, a preferential 
expansion of lighter single stars in the vicinity of the bi- 
nary. (Energy conservation is roughly applicable because 
the mass segregation timescale is shorter than the local 
relaxation timescale, by a factor of M/m, where M is 
the mass of the segregating object and m is the mass of 
a back ground star.) This e ffect is not included in the 
code of Ilvanova et alj (|2005T) . and is likely an important 
factor in the discrepancy between their results and ours. 

Another important factor, as suggested by Figs. [5] 
and 21 is that the lon g lived, high d e nsity binary burn- 
ing phase assumed by Ilvanova et al.l (|2005f ) may not be 
generic for globular clusters. Instead, the "core contrac- 



tion" phase may last a Hubble time, and the cluster cores 
we o bserve now ma y have been much less dense in the 
past (Frcgcau 2008) . Although the central density in our 
models increases steadily with time, the local density at 
the half-mass radius decreases with time, resulting in fi- 
nal half-mass relaxation times that are a factor of ~ 3 
longer than their initial values. In the cases where our 
models do enter the binary burning phase before a Hub- 
ble time, we find that the core binary fraction in this 
phase steadily decreases with time. This behav ior is con- 
sistent with the results of Ilvanova et alj (|2005l ). 

While the core binary fraction in the majority of our 
models increases with time, the overall cluster binary 
fraction remains roughly c onstant with t ime. This is in 
good agreement with the iHurlev et al.l (|2008f ) iV-body 
models inspired by NGC 6397, and supports their use of 
the currently observed binary fraction near the half-mass 
radius as a measure of the primordial binary fraction (al- 
though the validity of comparison with NGC 6397 is not 
obvious, since the iV-body models end with a factor of 
5 to 10 fewer stars than NGC 6397 currently contains). 
For the low-density /& = 0.05 model just described, 39% 
of the initial binary population remains at 14 Gyr, 43% 
escape the cluster due to the tidal field (compared with 
the 60% of single stars that escape in the same fashion), 
9% are destroyed via strong dynamical interactions of 
binaries, and 8% are destroyed via binary stellar evolu- 
tionary processes (possibly perturbed by dynamics). In 
other words, in this case the overall binary fraction re- 
mains roughly constant with time due to a balance be- 
tween preferential tidal stripping of single stars in the 
outskirts and preferential destruction of binaries in the 
cluster core. 

3. HIDING BINARIES 

When using the offset main sequence method, what 
observers measure is in fact the number of MS-MS bina- 
ries with mass ratios q > 0.5 relative to the total number 
of objects appearing in the main sequence (which may 
include apparent single MS stars, comprised of a MS star 
plus dim compact object companion). This fraction is 
then corrected to account for the low mass ratio MS-MS 
binaries that blend into the single MS, by adopting an as- 
sumed mass ratio distribution. This final corrected figure 
is what is usually quoted as the "observed binary frac- 
tion." However, there is no a priori reason to believe this 
quantity reflects t he underlying binar y fraction among 
stars of all types. IHurlev et al.l (|2008f ) showed that, for 
the low binary fraction cluster models they considered 
(fb ;$ 10%), the observed binary fraction is a good mea- 
sure of the true binary fraction in the outer regions of a 
cluster, but can be a serious overestimate in the core. 

Since the general nature of the relationship between 
the observed and true binary fraction is not obvious, 
we have plotted in Fig. [H] the evolution of our mod- 
els in core number density-o&serued core binary frac- 
tion space. The observed binary fraction is calculated as 
^Vms-ms/O^ms-ms + Nms + N M s-co), where N M s-ms 
is the number of MS-MS binaries of any mass ratio in 
the core, Nms 1S the number of single MS stars in the 
core, and -/Vms-CO is the number of MS-compact object 
binaries in the core that appear near the MS. We count 
a MS-compact object binary as near the MS if the total 
luminosity of the binary is less than 10% more than that 
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Fig. 6. — Same as Figure[3] but for main sequence binaries. 

of the MS star (this corresponds to a magnitude increase 
of 0.1), and if the luminosity- weighted temperature of 
the binary is less than 10% different from that of the 
MS star. Like the true binary fraction plotted in Fig. [3l 
small initial binary fraction models evolve toward larger 
binary fractions. However, large initial binary fraction 
models evolve toward drastically smaller observed binary 
fractions. As a relatively extreme example, the model 
with an initial binary fraction of 75% and initial cen- 
tral density of ~ 10 pc~ 3 has an observed core binary 
fraction of just 33% at the end of the simulation. The 
true core binary fraction at the end of the simulation 
is 91%. Of the core binaries, 23% are MS-MS binaries, 
32% are compact object-compact object binaries, and 
44% are MS-compact object binaries (see Table [T] for 
more details). As expected, the discrepancy between the 
observed and the true binary fraction is due to compact 
object-compact object binaries not being counted in the 
observed tally, and MS-compact objects masquerading 
as single stars on the main sequence. 

4. THE IMPORTANCE OF SOFT BINARIES 

An initial population of binaries that contains a sub- 
stantial soft component can be a significant cluster en- 
ergy sink, since the soft binaries are destroyed in dynam- 
ical scattering interactions. The result is that the core of 
a cluster born with many soft binaries will quickly con- 
tract as those binaries are ionized. Could soft binaries 
increase the concentration of a cluster so much that it 
would become core collapsed? 

The total energy in soft binaries, for a distribution flat 
in the log of the semimajor axis, is simply 

£ _ Nb(Eb,i is — i?b,amax) ^ -W&-Eb,hs /jx 
ln(fl max /(iniin) hl(dmax/ Gmin) 

where a m i n and a max are the limits on the semimajor axis 
distribution, -Eb.hs is the energy of a binary at the hard- 
soft boundary, -&b,amax is the energy of the least-bound 
binary, and Nf, is the total number of binaries. Assuming 
for simplicity a cluster of equal-mass objects (binaries in 



TABLE 1 
Population breakdown of 
core binaries at 14 gyr for 

THE ~ 10 2 ' 5 pc~ 3 INITIAL CORE 
DENSITY, 75% INITIAL BINARY 
FRACTION MODEL. 



type 


number 


fraction 


MS-WD 


139 


44% 


WD- WD 


101 


32% 


MS-MS 


71 


23% 


NS-WD 


1 


0.3% 


HG-WD 


1 


0.3% 



Note. — The third column 
is the fraction of the total num- 
ber of core binaries represented 
by that binary type. "MS" 
denotes main sequence, "WD" 
denotes white dwarf, "NS" de- 
notes neutron star, and "HG" 
denotes Hcrtzsprung gap star. 

this case) of mass m avc with mean ID velocity dispersion 
<t, this becomes 

fA b m ave g 2 
^b,s ~ 7-7 1 r • (/J 

ln(^(2 max / ttmin ) 

From the virial theorem, the total mechanical energy of 
a cluster is simply 2£ c ius = — ^Nm^c® 2 , where N is the 
number of cluster objects, so 

|-Ecius| iVln(a max /a min ) 

For an admittedly optimistic binary fraction of 1 (Nf, = 
N), and realistic binary semimajor axis limits of a m in = 
5 x 10~ 2 AU (corresponding to a contact binary during 
the pre-main sequence phase) and a max = 10 3 AU (corre- 
sponding to a 10 7 day orbital period), the energy in soft 
binaries is r* 10% of the total cluster mechanical energy! 

The question, of course, is if this amount of energy is 
sufficient to make a cluster concentrated enough to ap- 
pear to be core collapsed. For our working definition 
of core collapse we assume that a cluster core can be 
resolved with HST if its radius is at least 1 arcsecond 
in size. At a typical cluster distance of 10 kpc, this 
corresponds to ~ 0.05 pc. Starting with a King model 
of a given mass, binary fraction, central concentration 
Wo, and half-mass radius 77, , we calculate the total me- 
chanical energy of the cluster with in the half-mass ra- 
dius, Eh (|Binnev fc Tremainell2008[ ). We then calculate 
the energy of the soft binaries, £"b,s- This energy will 
be absorbed from the cluster when those binaries are 
destroyed in dynamical interactions in and around the 
cluster core. Keeping 77, fixed (since the timescalc for 
destruction of soft binaries is shorter than the half-mass 
relaxation time), we then calculate a new King model 
with half-mass energy E' h = E^—E^a (note that E'b.s > 
by construction). For the new King model we calculate 
the new central velocity dispersion and hence the new 
hard-soft boundary (which has moved to smaller binary 
semimajor axis), calculate the energy available in the 
newly soft binaries, and iterate until the solution con- 
verges. For a 5 x 10 5 M Q cluster with half-mass radius 
77j = 5pc and initial core radius r c = 1.9 pc (Wo = 6, 
concentration c = log 10 (r t /r c ) = 1.25), an initial binary 
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Fig. 7. — Evolution of the core radius and core binary fraction for 
our "high density" initial model with a 90% initial binary fraction, 
including soft binaries. The core contracts rapidly at the start of 
the evolution due to the destruction of soft binaries, and quickly 
enters the binary burning phase. 



fraction of 100% with semimajor axis distributed flat in 
log a from 5 x 10~ 2 AU to 10 3 AU is sufficient to drive the 
cluster to a Wo = 10, c = 2.3 King model with core ra- 
dius r c — 0.16 pc. (A Wo = 10 King model has maximal 
binding energy within Th for fixed rh and mass.) This is 
quite close to core collapsed, and may even be classified 
as such if viewed with a ground-based telescope. In fact, 
iWivanto et alj (|1985f ) showed that clusters enter the self- 
similar stage of evolution (the "onset" of core collapse) 
when Wo > 7.4, so such a model would reach core col- 
lapse quickly. We have repeated this calculation with a 
binar y distribution that is l o g-nor mal in orbital period, 
as in iDuquennov fc Mavoil (|199lD or iFischer fc Marcvl 
(|1992fl . with (logioPd) = 4.8 and <7i ogloPd = 2.3, where 
Pd is the period in days, and with the same limits on 
semi-major axis as above. The results are unchanged 
with this binary distribution, largely because its peak lies 
at wider orbits than the hard-soft boundary for globular 
clusters. 

To test this scenario numerically, we have run mod- 
els with a binary distribution extending well beyond the 
hard-soft boundary, to P = 10 7 d. Our "high density" 
model (cluster mass 9 x 10 4 M Q , standard wide mass 
spectrum, ry t = 1.0 pc initially) with /j = 0.9 (including 
soft binaries), evolves from /b.c = 0.9 and r c = 0.6 pc to 
/b, c = 0.4 and r c = 0.05 pc in just 3 Myr (see Fig. [7J. 
This is in striking agreement with the energy argument 
above, which predicts rapid evolution to r c = 0.1 pc for 
this model. Note that the energy argument assumes all 
soft binaries will be destroyed on a short timescale. To 
achieve this in practice requires efficient mass segregation 
of binaries into the core, which has been aided in this case 
by a wide mass spectrum, at the expense of inaccuracy 
in calculating i?b,s- After the rapid initial contraction of 
the core, the cluster quickly (after a few Myr) enters into 
a long-lived binary burning phase. 

From the preceding discussion, it is evident that the 
dynamical importance of soft binaries should not be ig- 
nored. If a cluster is born with significant numbers of soft 
binaries, its evolution may be vastly different from a sim- 
ilar cluster containing only hard binaries. First, the rate 
of binary destruction is greatly enhanced in clusters con- 
taining soft binaries, yielding a starkly decreasing binary 
fraction with time. Second, the binary burning phase 
is reached quickly (within a few Myr) due to soft binary 



destruction. When only hard binaries are present, the bi- 
nary burning phase may not be reached within a Hubble 
time, as shown for example by Fig. [51 The implications 
for our understanding of the current dynamical states of 
Galactic globular clusters are profound, as certain prop- 
erties of clusters can be explained by the majority of 
clusters currently being in the initial "core contrac tion" 
phase, and not yet in binary burning (|Fregeaull2008D . We 
have provided here just a cursory analysis of the effects 
of soft binaries. A more detailed study should certainly 
be undertaken in the future. 

5. DISCUSSION 

Independent of the details, it seems clear that the hard 
binary fraction in the core of a dense stellar system will 
generally increase with time (with the exception of an 
initial hard binary fraction > 90%). Yet there is no 
compelling evidence that clusters should be born with 
binary fractions smaller than the typical field value of 
~ 50%, and observations yield core binary fractions of 
just ~ 10% in Galactic globular clusters. If the observa- 
tions are to be taken at face value, how then can they 
be consistent with large init ial binary fraction s? One 
possibility, as pointed out bv lDavis et all (|2008D . is that 
the binary f raction is a strong function of primary mass 
()Ladal[2006l ). with the single star fraction increasing to 
~ 75% for M dwarfs and lighter stars. A iKroupa et al.l 
(| 19931 ) IMF with a 25% binary fraction from 0.1 to 0.5 
Mq and a 50% binary fraction from 0.5 to 100 M© yields 
an overall binary fraction of just 32%. 

Another possibility is that most binaries born in clus- 
ters are soft relative the cluster velocity dispersion, in 
which case they will be destroyed very quickly by dy- 
namics. If the binary period distribution is unifor m in 
logP from 0.1 to 10 7 d as in llvanova et al.l (|2005h . the 
32% overall binary fraction just suggested corresponds 
to a hard binary fraction of merely ~ 10% for a cluster 
with central density 10 6 pc -3 . As demonstrated above, 
the early, rapid destruction of soft binaries may lead to 
a binary-burning phase within a short time (< 5 Myr, 
depending on initial conditions). 

Aside from the initial binary properties, could it be 
that observations are under-counting the binary fraction 
significantly? When we measure the binary fraction us- 
ing an offset main sequence method similar to what ob- 
servers use, we find that clusters with large initial binary 
fractions (/& > 0.5) evolve toward smaller observed core 
binary fractions (/& < 0.5). The discrepancy between the 
observed and true core binary fractions is caused by com- 
pact object-compact object binaries not being counted in 
the sample, and MS-compact object binaries masquerad- 
ing as single stars. 

Could a binary be sufficiently wide to be resolved as 
two single stars and hence missed as a binary? For the 
wide-field camera on HST, one requires two turnoff mass 
stars in a binary to be separated by roughly 4 pixels for 
the binary to be resolved. For a cluster at a distance of 
10 kpc, this corresponds to a binary separation of ~ 4 x 
10 3 AU. For a cluster with velocity dispersion 10 kms -1 , 
this corresponds to a binary hardness of Gm/av^ « 2 x 
10 -3 , which is too soft to survive dynamically for even a 
short time. 

6. SUMMARY 
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Wc have described in detail our inclusion of the BSE 
single and binary stellar evolution ro utines in our Monte 
Carlo globular cluster evolution code (jHurlev et al.ll2000l . 
l2002f) . We have compared with the results of direct N- 
body simulations and found good agreement, suggesting 
that our implementation of BSE in our code is consistent 
with that in iV-body. 

We have used our newly upgrade Monte Carlo code 
to study the evolution of the core hard binary frac- 
tion in star clusters, and in particular attempt to set- 
tle the apparent disagreement between direct .ZV-body 
and simplified Monte Carlo techniques on its evolu- 
tion. Wc find that the core binary fraction generally 
increases with time, even for low initial core density mod- 
els (n c w 10 2 ' 5 pc -3 ), with only very small initial binary 
fraction models (/& < 0.05) producing the core binary 
fractions of ~ 10% observed today. The increase in the 
core binary fraction with time can be understood as an 
imbalance between mass segregation of binaries into the 
core (and single stars out of the core) and the destruc- 
tion of binaries in the core directly via strong dynamical 
encounters, and indirectly via dynamical perturbation of 
binary stellar evolution processes. The overall cluster 
binary fraction remains roughly steady with time, due 
to the additional effect of preferential tidal stripping of 
single stars from the cluster outskirts. 

This evolution, however, refers to the true binary frac- 
tion. When measuring the core binary fraction using 
an offset main-sequence method analogous to what ob- 
servers use, we find that the observed core binary fraction 
can seriously underestimate the true core binary fraction. 
This results from compact object-compact object bina- 



ries not being counted and MS-compact object binaries 
masquerading as single stars in the observed tally. In the 
course of creating more detailed models of M 67, 47 Tuc, 
M 4, and NGC 6397 to be compared with observations, 
we are now developing a data reduction pipeline that in- 
cludes simulations of spectra for every star. Among our 
near future plans is the creation of a cluster sky map in 
different bands, to which we can apply the MS binary de- 
tection method to determine more accurately how many 
binaries are missed by the method. 

Most of our discussion concerned hard binaries. How- 
ever, wc also considered the effects of a substantial pop- 
ulation of soft binaries. We found that the energy sink 
represented by soft binaries (for a typical binary period 
distribution) is sufficient to cause the core of a typical 
globular cluster to contract significantly. The result is 
not only a rapid, efficient destruction of a significant 
number of binaries at early times, but also a much earlier 
onset of the binary burning phase, resulting in enhanced 
binary destruction in the core with time. 
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